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1 Introduction 

In this paper, we compute reachable sets of differential inclusions, 

x(t) G F(x(t)), x(0) = x , (1) 

where F is a continuous set-valued map with compact and convex values. A solution of 
the differential inclusion is an absolutely continuous function x : [0, T] — > W 1 , such 
that for almost all t G [0,T], x(-) is differentiable at t and x(t) G F(x(t)). The solution 
set S T (x ) C C([0,T],R n ) is defined as 

S T {x ) = {x(-) G C([0,T],M n ) | x(-) is a solution of x(t) G F(x(t)) with x(0) = x }. 

The reachable set at time t, R(xo,t) C W 1 , is defined as 

R(x ,t) = {x(t) ER n \x(-) G SiCxo)}- 

In particular, we are interested in higher-order method for computation of a rigorous 
over-approximation of the reachable set of a differential inclusion. 

Differential inclusions are generalization of differential equations having multivalued 
right-hand sides, see [2], [8], [25]. They give a mathematical setting for studying differ- 
ential equations with discontinuous right-hand sides. In fact, taking a closed, convex hull 
of the right-hand side, one obtains a differential inclusion. Solutions of this differential 
inclusion are known as Fillipov solutions of the original differential equation; see [12] . 

One important application area for differential inclusions is control theory. Suppose 
we are given an interval [0, T], and absolutely continuous function x(-) which satisfies the 
inclusion ([T]), where F(x) = f(x,U) = {J u( zu f(x(t),u) for almost all t G [0, T\. It is 
known that if the set U is compact and separable, / is continuous, and f(x, U) is convex 
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for all x, then there exists a bounded measurable function u(t) G U, known as admissible 
control input, such that x(t) is the solution of the control system, 

x(t) = f{x{t),u{t)), x{0)=x . (2) 

The proof of the above is given in [2] , and with slight changes in the assumptions, also in 
[22] or [IS]. On the other hand, it is easy to see that each solution of a control system ^ for 
a given admissible control input is also a solution of a differential inclusion Therefore, 
if a control system is not completely controllable one may want to compute reachable 
sets corresponding to all possible inputs (u(t) 6 U) which is equivalent to computing a 
reachable set of a differential inclusion. 

Similarly, we obtain a differential inclusion from a noisy system of differential equations 

x{t) = f{x{t),v{t)), x{0)=x Ol v{t)eV. (3) 

Although the form of ^ and ^ are identical, the interpretation is different; in ([2]), the 
input u(t) can be chosen by the designer, whereas in ([3]), the input is determined by the 
environment. 

Differential inclusions can also arise as reduced models of high- dimensional systems 
of differential equations. For example, suppose we have a large-scale system given in the 
form of differential equation x(t) = f(x(t)). In general, it is very hard to analyse large- 
scale systems and most of the times performing model reduction is necessary. This gives 
a simplified model in the form of z(t) = h(z(t)) + e(i), where |e(£)| < e represents the 
error that occurred while simplifying the model. 

For reliability purposes many engineering systems require availability of verification 
tools. In order to verify a system, we must guarantee that an approximate solution will 
contain the actual solution of the system. If there is uncertainty in the system, lack 
of controllability, or just a variety of available dynamics, one needs to use differential 
inclusion models. For verification purposes, one needs to compute over-approximations 
to the set of solutions. 

An important tool in the study of input-affine control systems ^ is based on the Fliess 
expansion [T3], in which the evolution over a time-step h is expanded as a power-series 
in integrals of the input. A numerical method based on this approach was given in [15] . 
The method cannot be directly applied to study noisy systems ([3]), since for this problem 
we need to compute the evolution over all possible inputs, and this point is only briefly 
addressed. 

The first result on the computation of the solution set of a differential inclusion 
was given in [23], who considered Lipschitz differential inclusions, and gave a polyhe- 
dral method for obtaining an approximation of the solution set S(xq) to an arbitrary 
known accuracy. In the case where F is only upper-semi continuous with compact, convex 
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values, it is possible to compute arbitrarily accurate over- approximations to the solution 
set, as shown in [5]. 

Some different techniques and various types of numerical methods have been proposed 
as approximations to the solution set of a differential inclusion. For example, ellipsoidal 
calculus was used in [20], a Lohner-type algorithm in [TF], grid-based methods in [23] 
and [1], optimal control in [3] and discrete approximations in 0, HD1 ELI] , [II]- However, 
these algorithms either do not give rigorous over-approximations, or are approximations 
of low-order (e.g. Euler approximations with a first-order single-step truncation error). 
Essentially, the only algorithms mentioned above that could give arbitrary accurate error 
estimates are the ones that use grids. However, higher order discretization of a state 
space greatly affects efficiency of the algorithm. It was noted in [I] that if one is trying to 
obtain higher order error estimates on the solution set of differential inclusions then grid 
methods should be avoided. 

In order to provide an over-approximation of the reachable set of ([I]), we compute 
solutions of an "approximate" system 

y(t) = f(y(t),w k (t)), y(t k ) = x(t k ), w k {-) e W, 

for t G [t k ,t k+ i], and add the uniform error bound on the difference of the two solutions. 
We provide formulas for the local error based on Lipschitz constants and bounds on 
higher-order derivatives. The method is based on a Fliess-like expansion, and extends the 
results of [15] by providing error estimates which are valid for all possible inputs. 

We can obtain improved estimates by the use of the logarithmic norm. The logarithmic 
norm was introduced independently in [7], and [19] in order to derive error estimates to 
initial value problems, see also [26J. Using the logarithmic norm is advantageous over the 
use of Lipschitz constant in the sense that the logarithmic norm can have negative values, 
and thus, one can distinguish between forward and reverse time integration, and between 
stable and unstable systems. The definition of the logarithmic norm and a theorem on 
the logarithmic norm estimate is given in Section [2] 

The numerical result given in Section [6] were obtained using the function calculus 
implemented in the tool Ariadne [1] for reachability analysis and verification of hybrid 
systems. In particular, we use polynomial models for the rigorous approximation of con- 
tinuous functions. Polynomial model expresses approximations to a function in the form 
of a polynomial (defined over a suitably small domain) plus an interval remainder, and 
are essentially the same as the Taylor models of |24|. 

The paper is organized as follows. In Section [2j we give key ingredients of the theory 
used. In Section [3j we give mathematical setting for obtaining over-approximations of 
the reachable sets of a differential inclusion, and propose an algorithm. In Section |4| we 
consider differential inclusions in the form of input-affine systems. We derive the local 
error, give formulas for obtaining the error of second and third orders, and show how to 
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obtain the error of higher-orders. We extend the idea of obtaining over-approximations 
for input-affine systems to more general differential inclusions in Section [5] A numerical 
example is given in Section [6] We conclude the paper with a discussion on the theory 
proposed in Section [7j 



2 Preliminaries 

Below we give several results on differential inclusions and the computability of their 
solutions. For further work on the theory of differential inclusions see [2J, [8], [23], for 
computability theory see [28], and for results on computability of differential inclusions 
see [23], [5]. 

We canonically use the supremum norm for the vector norm in M. n , i.e., for x G M. n , 
\\x\\oo = max{|xi|, \x n \}. The corresponding norm for functions / : D C M. n — >■ R is 
sup x . g £, [| ^(ac) (| oo- The corresponding matrix norm is 
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n 

HQlloo = max { Y] \q ki \ \. 

k=l,...,n I.*—' J 



Given a square matrix Q and a matrix norm || • ||, the corresponding logarithmic norm 

is 

X(Q) = to + -' . 

There are explicit formulas for the logarithmic norm for several matrix norms, see [T6] . 
[7] . The formula for the logarithmic norm corresponding to the uniform matrix norm that 
we use is 

Aoo(Q) = m&x{q kk + \q ki \}. 

k ' J 

The following theorem on existence of solutions of differential inclusions and its proof 
can be found in [8]. Also, a version of the theorem and its proof can be found in [2J. 

Theorem 1. Let Del" and F : [0,T] x D =$ M. n be an upper semicontinuous set- 
valued mapping, with non-empty, compact and convex values. Assume that \\F(t, x))\\ < 
c(l + ||x||) ; for some constant c, is satisfied on [0,T]. Then for every xq G D, there exists an 
absolutely continuous function x : [0,T] — > M. n , such that x(to) = xq and x(t) G F(t,x(t)) 
for almost all t G [0, T]. 

A result on upper-semicomputability of differential inclusions was presented in [3]. 

Theorem 2. Let F be an upper- semicontinuous multivalued function with compact and 
convex values. Consider the initial value problem x G F(x), x(0) = Xq, where F is 
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defined on some open domain V C K n . Then the solution operator xq h- > St(xo) is upper- 
semicomputable in the following sense: Given an enumerator of all tuples (L, M\, M m ) 
such that F(L) C U™ 1 Mj ; it is possible to enumerate all tuples (I, J, K%, K k ) where 
I, K\, ...,K m are open rational boxes and J is an open rational interval such that for every 
x G I, every solution £ with £(0) = xo satisfies £(J) C U* =1 iiQ. 

In other words, it is possible to approximate the reachable sets arbitrarily accurately 
given a description of the differential inclusion and an arbitrarily accurate description of 
the initial state. 

The basic construction of our algorithm is based on the following theorem. The theo- 
rem and the proof can be found in [21 Corollary 1.14.1]. 

Theorem 3. Let f : X x U — >• X be continuous where U is a compact separable metric 
space and assume that there exists an interval I and an absolutely continuous x : I — >• R n ; 
such that for almost all t G I, 

x(t) G f(x(t),U). 

Then there exists a Lebesgue measurable u : I — » U such that for almost all t G I, 

x{t) = f{x{t),u{t)). 

We shall need the multidimensional mean value theorem, which can be found in stan- 
dard textbooks on real analysis, e.g., see [2Zj- We use the following form of the theorem. 

Theorem 4. Let V CM. 71 be open, and suppose that f : IR n — > W 71 is differ entiable on V. 
If x,x + h and L(x; x + h) C V , i.e., line between x and x + h belongs to V , 

f(x + h)-f{x)= [ Df{z{s))ds -h 
Jo 

where Df denotes Jacobian matrix of f , z(s) = x + sh, and integration is understood 
component-wise. 

The following theorem on the logarithmic norm estimate is taken from [16J. 

Theorem 5. Let x(t) satisfy differential equation x(t) = f(t,x(t)) with x(to) = Xq, 
where f is Lipschitz continuous. Suppose that there exist functions l(t), 5(t) and p such 
that X(Df(t,z(t))) < l(t) for all z{t) G cony{x(t),y(t)} and \\y(t) - f(t,y(t))\\ < 5(t), 
H^o) — 2/(^o)|| < P- Then for t > t we have 

\\y(t) -x(t)\\ < e f to l{s)ds (f> + J t e- f *o l ^ dr 5(s)ds S j . 
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In order to numerically compute the reachable set of a differential inclusion, we need a 
rigorous way of computing with sets and functions in Euclidean space. A suitable calculus 
is given by the Taylor models defined in [2Tj : 

Definition 6. Let / : D C M? — > K be a function that is (n + 1) times continuously 
partially differentiable on an open set containing the domain D. Let Xq be a point in D 
and P the n-th order Taylor polynomial of / around xq. Let I be an interval such that 
Then (p, I)p is called a Taylor model for / if 

f(x) — p(x — Xq) G J for all x G D 
Then we call the pair (P, /) an n-th order Taylor model of / around Xq on D. 

In Ariadne, we allow arbitrary polynomial approximations, and not just those defined 
by the Taylor series. We take p to be a polynomial on the unit domain [— and 
pre-compose p by the inverse of the affine scaling function s : [—1,+!]" — > D with 
Si(zi) = TiZi + mi. Instead of using an interval bound for the difference between / and p, 
we take a positive error bound e. We say (s,p,e) is a scaled polynomial model for / on 
the box domain D if s : [—1, +l] v — > D is an affine bijection and 

sup \ f(x) — < e. 

xeD 

In the special case D = [— l,+l] v , the unit box, we speak of a unit polynomial model 
(p, e) satsfying sup z6 r 1)+1 i 1) \f(z) — p(z)\ < e. We use the notation p o s -1 ± e to denote 
the polynomial model (s,p,e). 

Polynomial models support a complete function calculus, including the usual arith- 
metical operations, algebraic and transcendental functions. Formally, if op is an operator 
on functions, then there is a corresponding operator op on polynomial models satisfying 
the property that if /j are polynomial models for /j, i — 1, . . . ,n on common domain D, 
then op(/i, . . . , f n ) is a polynomial model for op(/i, . . . , f n ) on D. A full description of 
polynomial models as used in Ariadne is given in [B]. 

For the calculations described in this paper, it is sufficient to consider sets of the 
form S = f(D) for D = [— l,+l] m and / : lR m -» lR n . If pi ± ti are unit polynomial 
models for fi, then 

Sc S = p([-l,+l} m )±e 

= {x <E R n \ Xi = pi(z) + di for some z G [-1, +l] m and d G M n , \di\ < e;}. 

Here, p : [— l,+l] m — > 1R™ is the polynomial with components pi, and ±e is the set 
riiiit - e ii + e i]- The set S is an over- approximation to S. Note that by defining polyno- 
mials qi(z, w) = Pi(z) + eiWi, we have 

S = p([-l, +l] m ) ± e C +l] m+n ) 

yielding an over-approximation as the polynomial image of the unit box without error 
terms. 
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3 Approximation Scheme 

We consider differential inclusions in the form of noisy differential equations 

x(t) = f(x(t),v(t)), v(t)eV, (4) 

where x : K — > M. n , v(-) is a bounded measurable function, V C lR m is a compact convex 
set, / is continuous and f(x, V) is convex for all x G M. n . In order to compute an over- 
approximation to the reachable set of Q, we compute solution set of a different (an 
approximate) differential equation and add the uniform error bound on the difference of 
the two solutions. 

3.1 Single-step approximation 

Given an initial set of points Xq, define 

R(X , t) = {x(t) | x(-) is a solution of Q with x(0) G X } (5) 

as the reachable set at time t. 

Let [0, T) be an interval of existence of (JiJ). Let = t Q , ti, . . . , t n _i,t n = T be a 
partition of [0,T], and let h k = t k+l -t k . For x eR n and v(-) G L°°([t k , t k+1 ); R m ), define 
4>(x k ,v(-)) to be the point x(t k+1 ) which is the value at time t k+ i of the solution of ^ 
with x{t k ) = x k . 

At each time step we want to compute an over-approximation -R^+i to the set 

ieach(R k ,t k ,t k+1 ) = {(f)(x k ,v(-)) \ x k e R k and v(-) G L°°([t k , t k+1 ]; R m )}. 

Since the space of bounded measurable functions is infinite-dimensional, we aim to 
approximate the set of all solutions by restricting the disturbances to a finite-dimensional 
space. Consider a set of approximating functions W k C C([t k ,t k+ i];M. m ) parameterized 
as W k = {w(a k ,-) | a k G A C W}, such as w(a k ,t) = a ok + a lk (t - t k+1 / 2 )/h k where 
tk+1/2 —t k + h k /2 = {t k + t k+ i)/2. We then need to find an error bound e such that 

Vffc G L°°([t k ,t k+ i\; V), 3a k eAs.t. \\<f)(x k ,v k (-)) - (f)(x k ,w(a k , -))|| < e k . (6) 

Note that we do not need to find explicitly infinitely many a k s. Instead we need to choose 
the correct dimension (M. p ) and provide bounds on them to get desired error e k . Setting 
(p(x k ,a k ) = 4>(x k ,w(a k , ■)), i.e., (f> also denotes the solution of x{t) = f(x k ,w(a k , •)), with 
x(t k ) = x k , at t = t k+ i, we obtain the over-approximation 

R k+ i = {4>(x k ,a k ) + [-e k ,e k ) n \ x k G R k and a k G A}. 
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Define the approximate system at time step k by 

y{t) = f(y(t),w k (a k ,t)), y k = y(t k ), t G [t k ,t k+1 ]. (7) 

We would like to choose "approximating" functions w k = w(a k ,-) : [t k ,t k+ i] — > R, de- 
pending on x(t k ) and v (•), such that the solution of ^ is an approximation of high order 
to the solution of Q. The desired local error for this paper is at least of 0(h 3 ). Then we 
can expect the global error (cumulative error for the time of computation, [0,T]) to be 
roughly of 0{h 2 ). 

Without loss of generality, we assume that x{t k ) = y{t k ) for all k > 0. To be precise, 
initially, we assume x(t ) = y(to). After obtaining an over-approximation Ri, to the 
solution set at time t 1 , we use Ri at the set of initial points of both the original system (|4]) 
and its approximation ^ for the next time step. Thus we have x(ti) = y(ti) G R\. We 
compute i?2, and consider it to be the set of initial points for both equations at time t 2 . 
Proceeding like this, we have x(t k ) = y(t k ), for all k > 0. 

The local error for a time-step consists of two parts. The first part is the an- 
alytical error given by ([6]). The second part is the numerical error which is an in- 
terval remainder of the polynomial model (see Definition [6J representing the solution 
(p(x k ,a k ) of x{t) = f(x(t),w k (a k ,t)). We represent the time-t^ reachable set R k = 
{h k (s) + [— e k ,e k ] n \s G [— 1, +l] Pfc }, as a polynomial model whose remainder consists 
of both numerical and analytical error. Here, p k is the number of parameters used in 
the description of R k . The inclusion R(X ,t k ) C R k is guaranteed by this approximation 
scheme. 

Note that our method only guarantees a local error of high order at the sequence of 
rational points {t k } which is a priori chosen. If one is trying to estimate the error at 
times t k < t < t k+ i for any k along a particular solution, a different formula should be 
used such as a logarithmic norm estimate based on Theorem [5] 

3.2 Algorithm for Computing the Reachable Set 

In this section we present an algorithm for computation of the solution set of Q, using 
the single step computation presented earlier. 

Algorithm 7. Let R k = {h k (s) ± e& | s G [—1, +l] Pfc } be an over-approximation of the set 
R(X ,t k ). To compute an over-approximation of R(X ,t k+1 ): 

1. Compute the flow 4> k (x k ,a k ) of 

x(t) = f(x(t),w k (a k ,t)), x(t k ) = x k , 

for t G [t k ,t k+ i], x k G R k , and a k G A. 
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2. Compute the uniform error bound Ek for the error of approximating x = f(x(t),v(t)) 
by x = f(x(t),w(t)). 



3. Compute the set Rk+i which over-approximates R(xo,tk+i) as Rk+i 3 {4>(xk,ak) ± 
e fc | x fc G -Rfc, a fc G A}. 

4. Reduce the number of parameters (if necessary). 

5. Split the new obtained domain (if necessary). 

Step [l] of the algorithm produces an approximated flow in the form faix^, a^) m 
(f>(xk,w(ak, •)) which is guaranteed to be valid for all Xk G -Rfc- In practice, we cannot 
represent cf) exactly, and instead use polynomial model approximation with guaranteed 
error bound 0. In Step |ij we add the uniform error bound St to make sure an over- 
approximation is achieved. In Step [3j we compute a new approximating set by applying 
the approximated flow to the initial set of points to obtain a solution set Rk+i = {0(/i(sfc)± 
efc, afe) ±£fc}- Steps [5] and [5] are crucial for the efficiency and the accuracy of the algorithm, 
as explained below. 

It is important to notice that the number of parameters (a^ initially) grows over the 
time steps. At each time-step, the number of parameters doubles, unless certain reduction 
of parameters is applied. The easiest way to reduce the number of parameters is to replace 
the parameter dependency by a uniform error, but this can have a negative impact on 
the accuracy. Another way to reduce number of parameters is using orthogonalization, 
though this is only possible for affine approximations using currently known methods. 

It is also of importance to realize that if the approximating set becomes too large, it 
may be hard to compute "good" approximations to the flow and/or the error. In this 
case, we can split the set into smaller pieces, and evolve each piece separately. This can 
improve the error, but is of exponential complexity in the state-space dimension. 

4 Input- Affine Systems 

In this section, we restrict attention to the input-afnne system 

x(t) = f(x(t)) + ^2g i (x(t))v i (t); x(t ) = x . (8) 

For some r > 1 which depends on the desired order, we assume that 

• / : R n -> R n is C r function, 

• each gi : IR n — > W 1 is C r function, 
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• Vi(") is a measurable function such that Vi(t) G [—Vi, +Vi] for some Vi > 0. 
Then the equation ^ becomes 

m 

y(t) = f(y(t)) +^2gi(y(t))wi(a k ,t); y{t k ) = y k , t G [t k ,t k+1 ]. (9) 
i=i 

In what follows, we assume that we have a bound B on the solutions of (J8j) and (|9]) 
for all t G [0, T]. We take constants Vi, K, K^, L, Li, H, A such that 

< V t , \\f(z(t))\\ < K, \\ gi {z{t))\\ < Ki X(Df(-)) < A, 
\\Df{z{t))\\ < L, \\D gi (z(t))\\ < ^, ||£> 2 /(^(*))ll < H, ||D 2 ^(^(t))|| < H it 

for each i — 1, m, and for all t G [0, T], and z(-) G B. We also set 

m mm 

K' = Y,VKi, L' = Yy-iLi H^J^ViHi. 

1=1 1=1 1=1 

Here, Df denotes the Jacobian matrix, D 2 f denotes the Hessian matrix, and A(-) denotes 
the logarithmic norm of a matrix defined in Section [2j 

We proceed to derive higher order estimates on the error by considering several dif- 
ferent cases. In each of the cases, u>j(a, ■) is a real valued finitely-parametrised function 
with a G A C K . In general, the number of parameters N depends on the number of 
inputs and the order of error desired. 

In what follows, we write h k = t k+ \ — t k , t k+ i/ 2 = t k + h k /2 = (t k + t k+ i)/2, and 

?(*) = ft„ v( s ) ds - 

4.1 Error derivation 

The single-step error in the difference between x k+ i and y k+ \ is derived as follows. Writ- 
ing (|sj) and ([9]) as integral equations, we obtain: 

rtk+x ™ 

x{t k+l ) = x{t k ) + f(x(t)) + ^ gi (x(t))vi(t)dt; (11a) 

Jtk j = i 

y(t k+l )=y(t k )+ / f(y(t)) + J29i(y(t))wi(t)dt. (lib) 
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Since we can take x(tk) = y{tk) as explained in Section [3l we obtain 



x(t k+1 ) - y(t k+ i) 



-fc + 1 



f(x(t))-f(y(t))dt 



tk+l 



+ 9i(x(t))vi(t) - gi (y(t))wi(t)dt. 



(12a) 
(12b) 



Integrating by parts the term (12a), we obtain 



(12a) 



(t-t k+1/2 ){f(x(t))-f(y(t))) 



c fc+i 



tk+ \t-h +1/2 )j t {f(x(t))-f(y(t)))dt 



tk 



(h k /2)(f(x(t k+1 ))-f(y(t k+1 ))) 



tk + l 



\t -t k+1/2 )(Df(x(t))x(t) - Df(y(t))y(t))dt. 



There two ways that we deal with term (12b). First we rewrite the term inside the integral 

as 

gi (x(t))vi(t) - gi {y(t))wi{t) = (gMt)) - gMt))) Wi (t) + gi(x(t)) ( Vi {t) - Wi (t)), 
and then integrate by parts the second term to obtain 

S) = X) / - gMt))) Wi(t) dt 

1=1 

171 d 



r i f k+1 d ( 



i=l 
m ntk+i 

J2 / (9i(x(t))-gi(y(t)))wi(t)dt 
i=i 

m 

+ / J gi{x{tk+i)){vi{tk+i) - Wi(t k+1 )) 
i=i 

m rtk+i 

- / D gi (x(t)) x(t) (vi(t) - Wi(t)) dt 



gi{x{t))\ {vi{t)-Wi{t))dt 



(13a) 
(13b) 
(13c) 



i=i Jl k 

The second derivation is obtained just by integrating by parts, 
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( pbj ) = \9i(x(t))vi(t) - gMt))wi(t) 



i=l 



l k + l 



i=l " l * 



^ gi(x(t k+1 ))vi(t k+1 ) - gi(y(t k+1 ))wi(t k+ i) 



i=l 



m rt 



•fe+i 



D gi (x(t))vi(t)x(t) - D gi (y(t))wi(t)y(t) dt 



(14a) 
(14b) 



i=i " l * 



Equations (12a) and (13) can be used to derive second-order local error estimates. 



By applying the mean value theorem (Theorem [4]) we obtain 



f(x(t k+1 )) - f(y(t k+1 ) = I Df(z(s))ds (x(t k+1 ) - y(t k+1 )) 

'o 



Hence 



pa] ) = (h k /2) / Df(z(s))ds (x(t k+1 ) - y(t k+1 )) 



tk+l 



(t - t k+1/2 )(Df(x(t))x(t) - Df(y(t))y(t))dt. 



Separate the second part of the integrand in (15b) as 



Df(x(t)) x(t) - Df(y(t)) y(t) = Df(x(t)) (x(t) - y(t)) 

+ (Df(x(t)) - Df(y(t))) y(t) 

The first term of the right-hand-side can be expanded using 

rn 

x(t) - y(t) = f(x(t)) - f(y(t)) + ^2(gMt)) - g t (y(t)))w t (t) 

m 



i=l 



(15a) 
(15b) 



(16a) 
(16b) 



i=l 
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Hence we obtain 



pa) = (h k /2) / Df(z(s))ds (x(t k+1 ) - y(t k+1 )) 
Jo 

L+ \t ~ 4+i/ 2 ) Df(x(t)) (f(x(t)) - f(y(t))) dt 



m rt 



- Z) / (* - A/(z(t)) - gi{y{t)))wi{t) dt 

i=l 

m rtk+i 



i=l 



(* - 4+1/2) (Df(x(t)) - Df(y(t))) y(t)dt 



(17a) 
(17b) 

(17c) 

(17d) 
(17e) 



where (17a) is ( 15a[ ), (Jl7|3-d) come from (16a), and (17e) comes from (16b). Note that 
for any C 1 -function h(x) we can write 

\h(x(t)) - h(y(t))\ < \\Dh(z(t))\\ ■ \x(t) - y(t)\ 

where z(t) G conv{x(t), y(t)}. This will allow us to obtain third-order bounds for 



terms (|17p,c,e). In order to obtain a third-order estimate for term (17d), a further inte- 
gration by parts is needed. We obtain: 



8=1 



tk 



(B = - Z [Df(x(t)) g t (x(t)) [is- t k+l/2 )(v(s) - w(s))ds 

+1 (D 2 f(x(t)) 9i {x(t)) + Df{x{t))D gi {x{t))) x(t) 



;isd) 



{s - t k+1/2 ){vi{s) - Wi(s))ds dt. 



Using similar type of derivation as for the derivation of (17), again using the mean value 
theorem and integration by parts, we obtain 



Pa) + g4b|) = J2 D gi {z{s))ds (x(t k+1 ) - y(t k+1 ))wi(t k+1 ) 
i=i J ° 

m 

+ Z 9i( x k+l) (Vi(t k+ i) - Wi(t k+1 )) 
i=l 

m rtk+i 

-J2 (Dgi(x(t)) - D 9i {y{t))) y(t) w t (t)dt 



(19a) 
(19b) 
(19c) 



i=l 
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J2 / D 9l (x(t)) (f(x(t)) - f(y(t))) w t {t) dt (19d) 
i=i Jt * 
m rtk+i 

/ D 9i (x(t)) f(x(t)) (vi(t) - Wi(t)) (19e) 



i=l •'^ 



Dg i (x(t))(g j (x(t))-g j (y(t)))w j (t)w i (t)dt (19f) 

i=i j=i 

m m tk+1 

D ^{t)) 9Mt)) (Vjityiit) ~ WjtyWiit)) dt. (19g) 

i=l 7=1 



The term (19e) can be further integrated by parts to obtain 



§Uk$ = -J2[DgMt)) f(x(t)) J (v(s) - w(s))ds 



tk+1 

tk 



m rtk+i 

J2 / (D 2 gi(x(t)) f(x(t)) + D 9i (x(t)) Df(x(t)))x(t) fait) - Wi(t)) dt 
i=i Jt * 

(20e) 



and the term (19g) to obtain 



& = -J2J2\Dgi(x(t)) 9j (x(t)) [ ( 



Vj{s)vi(s) — Wj(s)wi(s))ds 



t=i j=i 



m m „i 



i=i j=i ^ (20g) 

(fj(s)'Oi(s) — Wj(s)wi(s))ds dt. 



Equations (17 20) can be used to derive third-order local error estimates. 



4.2 Local error estimates 

We proceed to give formulas for the local error having different assumptions on functions 
/(•), Qi(-) and Wi(-). We present necessary and sufficient conditions for obtaining local 
errors of 0(h), 0(h 2 ), 0(h 3 ), and give a methodology to obtaining even higher-order 
errors. In addition, we give formulas for the error calculation in several cases. 
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4.2.1 Local error of 0(h) 

Theorem 8. For any k > 0, and all i = 1, ...,m, if 

• /(•) is a Lipschitz continuous vector function, 

• gi(-) are continuous vector functions, and 

• Wi(t) = on [t k ,t k+1 ], 

then the local error is of 0(h). Moreover, a formula for the error bound is: 

p A/i fe _ i 



\x(t k +i) - y{h+i)\ < h k K' 



A hi 



(21) 



Alternatively, we can use 



\x(t k+1 ) - y(t k+1 )\ <h k [2K + K'). 



(22) 



Proof. Since Wi(t) = 0, we have y(t) = f(y(t)). Using the bounds given in (10), we can 
take l(t) = A in Theorem [Bl Further, and 



m mm 

y (t)-(f(y(t)) + Y,9i(y(t)Mt)) = 5>0/WMW) <Y,K i V i = K' 

i=l i=l i=l 



so we can take S(t) = YliLi^-i^i- Hence the formula (21) is obtained directly from 

Theorem [5) Note that (e Ahk - l)/(Ah k ) = 1 + Ah k /2 H is 0(1), so the local error is 

of 0(h). Equation ( |2~2~| ) can be obtained by noting that sup fg [ ffe 4fc+i ] \ \f(x(t)) — f(y(t))\ \ < 



4.2.2 Local error of 0(h 2 ) 

Theorem 9. For any k > 0, and all i — 1, m, if 

• /(•), gi(-) are C 1 vector functions, and 

• Wi(-) are bounded measurable functions defined on [t k ,t k+ i] which satisfy 



•fc+i 



Vi(t) — Wi(t) dt = 0, 



(23) 



then an error of 0(h 2 ) is obtained. 
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Proof. To show that the error is of 0(h 2 ), we use equations ( 12"fl3~ ). The equation (12a) 
is in the desired form, i.e., of 0(h 2 ), since we can write 



•fe+i 



f(x(t))-f(y(t))dt 



tk 



< 



h L sup te{tkjtk+l] \\x(t) -y(t)\\, 



and sup tg( - tfcjtft+i ) — y(t)\\ is of 0(h) by Theorem [5j Similarly, equations (13a) and (13c) 
are of 0(h 2 ). Note that the equation (13b) is zero due to (23). The theorem is proved. 

□ 

In order to be able to compute the errors, we need the bounds on the functions Wi(-). In 
particular, we can restrict Wi(-) to belong to certain class of functions, such as polynomial 
or step functions. 

Theorem 10. For any k > 0, and all i — 1, ...,m, if 

• /("); 9i(') OT C 1 vector functions, and 

• Wiit) are real valued, constant functions defined on [tk,tk+i] by Wi = f* k+1 Vi(t)dt, 
then a formula for calculation of the local error is given by 

A h k _ Y 



x(t k+1 ) - y(t k+1 )\\ <h 2 k ((K + K') L'/3 + 2 K' (L + ti) 



A hi 



(24) 



Proof. To derive (24), we obtain —y(tk+i)\\ from equations (12a) and (13). Using 

the bounds given in (10), it is immediate that ||x|| < K + Y^i+x an d straightforward 



to show that < V% and \vi(t) — Wi(t)\ < IVihk for t G [tk,tk+i]. However, we can 

get a slighly better bound \vi(t) — Wi(t)\ < Vih^/2 by considering the following: Without 
loss of generality, assume t G [0,h], and let 



a,i(t) 



Vi(s)ds, bi(t) 



h-t 



h-t 



vAs) ds 



and define 



Wi{t) = (tai(t) + (h-t)bi(t))/h. 



Then, Wi = Wi(t) is constant for all t G [0,h]. Notice that Vi(t) = ta(t) and Wi(t) 
(t/h)(ta(t) + (h- t)b(t)). Hence, we have 

Vi(t) - Wi(t) = t(h - t)(a(t) - b(t))/h, 
\viit) -Wi(t)\ =t(h-t)\a(t) -b(t)\/h < Vih/2. 
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Additionally, we can prove that J t k+1 \vi(t) - Wi(t)\dt < Vi h\/3. Take z(t) to satisfy the 
differential equation z(t) = f(z(t)). From Theorem |5j we have 



\x(t) - z(t)\\, \\y(t) - z(t)\\ < h k (TKM 



Ah k 



and hence 



D Ah k _ i 



Ah k 



\\x(t)-y(t)\\<2h k (j2 KiV \ 

4=1 

for t E [tfe, tfc+i]. Taking the norm of the equations ( 12a|13a|13c ) we obtain 

e Ahk - r 



x(t k+1 ) - y{t k+l ) || < [ L+1 (l + J2 ViL t ) (2 h k (V K i V i 



i- 4=1 
m 



4=1 



A/i* 



+ ^ ^ + V i K i) ~ I dt 



i=i 

m 



(l + E ^) 2 (E « ) + 5 (E (* + E »s«i) ■ 

4=1 v i=i fe / i=i j=i / 



Using if' and L', we get the desired formula (24). 



□ 



e An -l 



Remark 11. Note that as A — >• 0, then x h 
[5j In fact, if A = 0, we get 



— > 1. This is also consistent with Theorem 



and therefore, 



m 

\\x(t)-y(t)\\ <2h k (j2 K * V i) 

i=i 



\x(t k+1 ) - y(i fc+1 ) || < h\ ({K + K') L'/3 + 2 K' (L + L')) , 



(25) 



which is still of 0(h 2 ). Further, we will not give explicit formulas for the error when 
A = 0. 



Theorem 12. If all assumptions of Theorem 10 are satisfied, and in addition /(•) is C 2 , 
then a formula for calculation of the local error can be given by 

„Ah k 



1 - (h k L/2))\\x(t k+1 )-y(t k+1 )\\ < {h 2 k /3) \3K'L 



1 



Ahi 



+ L'{K + K') 



(h 3 k / A) K' (l U + L 2 + H(K + K') 
(ll/ifc/24) (H K' + LL')(K + K'). 



Ahk 
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Proof. The same bounds on functions apply as in Theorem 10 The formula for ||rr(ifc+i) — 
y(tk+i)\\ is then obtained by taking norms of terms in equations (17) and (13). □ 

Remark 13. The computation of the error bound is complicated by that fact that \vi(t) — 
Wi(t)\ is not uniformly small. This means that the terms g(x)(vi — Wi) must be inte- 
grated over a complete time step in order to be able to use the fact that f^ +1 Vi(t) dt = 

ft* +1 Wi(t) dt, and this must be done without first taking norms inside the integral. As 
a result, we cannot apply results on the logarithmic norm exactly directly. Instead, we 
"bootstrap" the procedure by applying a first-order estimate for \\x(t) — y(t)\\ valid for 
any t G [t k ,t k+1 ]. 

4.2.3 Local error 0(h 2 ) + 0(h 3 ) 

We can attempt to improve the error bounds by allowing Wi(t) to have two independent 
parameters. In the general case, we shall see that this gives rise to a local error estimate 
containing terms of 0(h 2 ) and 0(h 3 ), rather than the anticipated pure 0(h 3 ) error. 
We require Wi(t) to satisfy the equations 



Vi(t)-Wi(t)dt = Q; / (t-t k+1/2 ) ( Vi (t) - Wi (t))dt = Q. (26) 

Jt k 

If the Wi are taken to be affine functions, Wi(t) = a i + a i: i(t — t k +i/2)/h k , then we 
have 

I r-tk+i 12 /"**+! 

O'ifl — t~ I Vi(t)dt; a i,i = T2 / Vi(t)(t -t k+ i /2 )dt. (27) 
h k Jt, ri k J tk 



It is easy to see that 

Kol < ki| <3Vi, \wi(t)\ < 5Vi/2, and \w(t)\ < 3Vi/2h k (28) 
and it can further be shown that 

\a it i\<^V i {l-(a it0 /V i ) 2 ). (29) 
An alternative is to use step functions for u>j, such as 

di,0 if tk < t < tfc+l/2 



Wi(t) 



a 



;,i if t k+ i/2 < t < t k+ i. 



Then 



I rtk+i 4 rtk+i 

a ifl = TT / Vi(t)dt--^ I Vi(t)(t - tfc+i/a) dt 

hk Jt k h k J tk 

fli,i = T- / v i {t)dt+— 1 I Vi(t)(t - tjfe+1/2) dt. 

n k Jt k ri k Jt k 
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Hence 



(30) 



Theorem 14. For any k > 0, and all % — 1, ...,m, if 

• /(•) is C 2 vector function, 

• &(■) ° re non-constant C 2 functions, and 

• the Wi satisfy (|26|), 



then an error of 0{h 2 ) is obtained. Moreover, if the Wi are affine functions, W{(t) 
a i,o + a i,i(t — tk+1/2) /hk, then a formula for calculation of the error is given by 



l - L(h k /2) - h k L') \\x(t k+1 ) - y(t k+1 )\\ < (h 2 k /4)L' (11K + (69/2) AT') 



e Ah k _ 1 



+ (7/^/8) K' ((AH' + H)(K+ (5/2)K') + L 2 + ((9/2)L + 5L') L') 
+ (7^/48) (HK' + L L') (K + K') . 



Proof. With the assumptions of the theorem, we can improve the terms (17d) and (19e) 
such that they become (18d) and (20e), which are of 0(h 3 ). In addition to the bounds 
obtained in (28), we use 



x 



<K + ^K,V t = K + K' 



i=i 



){t)\\ <K + -Y J K t V i = K + (h/2)K l 



\x{t)-y{t)\\ < 



7h 



i=i 



e Ahk - 1 _ 7h k , e Ahk - 1 
Ah k ~ 2 Ah k 



The formula for the error, ||x(tfc+i) — y(tk+i)\\ with terms (18d) and (20e) is then easily 
obtained. The theorem is proved. □ 

We now show that with the assumptions of the theorem we cannot in general obtain 
an error of 0(h 3 ). Specifically, we assume that Wi(t) are two-parameter polynomial or 
step functions satisfying 



tk+i 



Vi(t) - Wi(t) dt 



tk+l 



(t-t k+1/2 )(v l (t)-w i (t))dt = 0. 



The following counterexample gives a system for which only 0(h 2 ) local error is possible. 
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Example 15. Consider the following input-affme system which satisfies assumptions in 
Theorem [14J 

±i = x 2 + v 1 + xiv 2 ; x 2 = x 1 + v 2 ; x(t k ) = x k . 

Take inputs 

vx(t) = sin - h) J , v 2 (t) = cos f^(* - *k)J • 



Using (26), we get wi(t) = — (6/tt hk)(t — ^+1/2), = 0. Therefore, an approximation 



equation looks like 

y 1 =y 2 + wi; y 2 = yi 
As shown in the previous section, the only term which might not have order h\ is the 



term in (19g) which is reduced to 

I W Dg 2 (x(t)) 9i (x(t)) Vi (t)v 2 (t)dt, 

1=1 

since Dgi = 0. When i = 2, we have |^(% 2 (^)) = Vi(t)vi(t), and hence we can integrate 
by parts once more to get the 0(h 3 ). Then we are left with 

Dg 2 {x{t))g 1 {x{t))v 1 {t)v 2 {t)dt = -^ [1 0] T , 



a term of 0(h 2 ). 



4.2.4 Local error of 0(h 3 ) 

We showed that for a general input-affine system, a local error of order 0(h 3 ) cannot be 
obtained using affine approximate inputs w(a, t). However, if in addition, we assume that 
gi(-) are constant functions or we have a single input then we can obtain a local error of 
0(h 3 ). If <7j(-) are constant functions, then the error calculation is equivalent to the error 
calculation of an even simpler case, so called additive noise case. The equation is then 
given by 

x(t) = f(x(t))+v(t). (31) 
Here, v(t) = (vi(t), v n (t)) is vector-valued. 

Corollary 16. For any k > 0, 

• if the system has additive noise, 
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/(•) is a C 2 function, and 



• Wi(t) are real valued functions defined on [t k ,t k +i] which satisfy equations (26), 

then an error of 0{h 3 ) is obtained. Moreover, for Wi(t) = a^o + a>i,i(t — t k+ \/2)/h k , the 
formula for the local error is given by: 



(1 - (h k /2)L) \\x(t k+1 ) - y(t k+1 ) || < ^ h 3 k K'H(K + K') 

+ - hi K' (L 2 + H(K + 5K'/2) 



n Ah k 



(32) 



The formula for the error in additive noise case is simplified because V — H' — 0. If we 



write ||f(t)|| = K', then the result follows directly from Theorem 14 
Corollary 17. For any k > 0, if 

• the input-affine system has single input, i.e., m — 1 in (J8 

• /(•) and g(-) are C 2 functions, and 



• w(t) is a real valued function defined on [tk,t k+ i] which satisfies equations (26), 

then an error of 0(h 3 ) is obtained. Moreover, for w(t) = a + a\(t — ^+1/2) , the formula 
for the local error is given by 



(1 - {h k /2)L - h k L') \\x(t k+1 ) -y(t 
7 hi 



k+l, 



< 



k K' ({H + 10 H'){K + (5/2)^) + L 2 + (25/2) L U + 25 (L') 2 ) 



A hi 



hi 



+ -±(K + K') ((7/Q)(H K' + LL') + 28 (H' K + LL') + 29 (H' K' + (L') 2 )) . 
48 

Proof. The result follows since the only term which is not 0(h 3 ) in ( 17|19 ) is ( |19g ). In 
the one-input case, this simplifies to 

■fe+i 

Dg{x{t)) g{x{t)) (v(t) v(t) - w(t) w(t)) dt. 

However, we can integrate by parts to obtain 

~Dg{x{t))g{x{t)){v{t) 2 -w{t) 2 ) 



(19g) 



tk 



tk+1 



D(Dg(x(t)) g(x(t))) x(t) (v(t) 2 - w(t) 2 ) dt. 



The first term vanishes since v(t kl ) = w(t k+ i), and the second is 0(h 3 ) since v(t) and w(t) 
are 0(h). Taking all the bounds as in Theorem 14, the formula is easily obtained. □ 
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Observing the error given by equations (17) and (19) , we see that if in addition to 



satisfying equations given in (26), the functions Wi(-) also satisfy 



■fc+i 



Vi(t)i!j(t) — Wi(t)u)j(t) dt = 0. 



(33) 



tk 



then we could get an error of 0(h 3 ). The question remains as to whether we can find 
functions lUj(') that satisfy the conditions ( 26~f3~3 ). Since the functions lOj(-) cannot be 
computed independently any more, the number of parameters of each Wi(-) will depend 
on the number of inputs. 

Theorem 18. For any k > 0, if 

• /(•)> 9i{~) are C 2 real vector functions, and 

• Wi(a,ifi, Oi jP -i, t) are real valued, defined on [tk,tk+i], and satisfy 



•fc+i 



Vi(t) - Wi(t)dt = 



*fc+i 



(t - t k+1/2 ) (Vi(t) - Wi(t)) dt = 



(34) 



■h+l 



Vi(t)iij(t) — Wi(t)vjj(t) dt = 0, 



for all i,j = 1, ...,m, then an error of 0(h 3 ) can be obtained. Note that it suffices to take 
j < i in (34), and that the number of parameters p in each W{ must satisfy p > (m + 3)/2. 
Taking polynomials of minimal degree d, we obtain d = \(m + 1)/2~|. 

Proof. If we can find Wi(t) that satisfies above, then it is obvious that the only remaining 
0(h 2 ) term (19g) can be integrated by parts once more in order to give a term of 0(h 3 ). 



This follows from Theorem Rl Corollary 17 and the formulae in Section 4.1 



To see that we can find the desired functions Wj(-), we consider polynomial approxi- 
mations Wi of degree d — p+1. We will show that it is possible to solve for the parameters 
of Wi's. If m - 
m + m + m(m 



1, see Corollary 17 The system of equations (34) consists of at most 
l)/2 = m(m + 3) /2 independent equations. To see that third equation 
in (34) has at most m(m — l)/2 independent equations necessary to be zero, notice that 
when i = j we have 



: fe+i 



Vi(t)vi(t) - Wi(t)wi(t) dt 



(l/2)[v 2 (t k+1 ) 



w 2 {t k+l )}, 



tk 
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# Inputs 


^Equations 


Degree 


^Parameters 


m 


m(m + 3)/2 


d 


m(d+ 1) 


1 


2 


1 


2 


2 


5 


2 


6 


3 


9 


2 


9 


4 


14 


3 


16 


5 


20 


3 


20 


6 


27 


4 


30 


10 


65 


6 


70 



Table 1: The number of independet equations which need to be solved, the minimal 
degree of a polynomial Wi(-) required, and number of available parameters in order to 
obtain 0(h 3 ) local error for m inputs. 

and therefore we can integrate by parts once more to get error of 0(h 3 ). When j > i 
integration by parts gives 

/ Vi(t)Vj(t) - WiftWjit) dt = [v^t) v 3 (t) - Wi^wtf)]** 1 

tk+1 

Vi(t)vj(t) — Wi(t)wj(t) dt 

and the first term vanishes since Vi(tk + i) = ■u)j(t/ c+1 ). The number of parameters that 
each Wi(-) has is p = d + 1. Thus, in total, we have mp parameters. In order to guarantee 
that we can solve all the equations for the tOj(-)'s, we need that mp > m(m + 3)/2. This 
implies that p > (m + 3)/2. Taking polynomials of minimal degree, we see that we require 
d=\(m + l)/2]. □ 

In what follows, we write C(n,m) = n\/{m\ (n — m)!), the formula for combinations 
(selecting m elements among n elements). 

In Table [TJ we present the degree of tOj(-) needed for one to obtain 0(h 3 ) for different 
number of inputs. In addition, the number of equations involved and the number of 
independent parameters in m functions that have to be found are given. 




4.2.5 Higher Order Local Error 

It is possible to generalize the approach used to generate 0(h 3 ) local error. With addi- 
tional smoothness requirements on the functions /(•) and <7i(-)' s > we can g e t even higher- 
order local errors. In order to simplify the notation, we set go = f and vq = 1. Then the 
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input-affine system (|8j) becomes 

m 

x(t) = ^2gi(x(t))vi(t). 
Let gi G C r for all i = 0, m, and denote by 

m 



i=0 



the corresponding approximate system. The local error of 0(h r+1 ) can be obtained if 
Wi(a,i,t) is finitely parametrised, = (a^o, a^) with d being sufficiently large, and 
satisfying 



-k+l 

Vi (t) dt 

tk+l ft 

Vj(t) / Vi(s) ds dt 

tk Jtk 
tk + l ft fS 

Vk(t) / v j( s ) / Vi(r) dr ds dt 
tk J t k Jt k 



-fc+i 



Wiit) dt 



(35a) 
(35b) 



tk+l ft 

Wj(t) / Wi(s)ds dt 
tk Jt k 

tk+l ft fs 

Wk{t) / Wj(s) / Wi(r) dr ds dt (35c) 

tk t k J t k 



tk + l fSr fS2 

V ir (s r ) / Vtr^Sr-x) ■ ■ V il (si)dSi • ■ ■ dSr^dSr = 

tk J t k Jtk 

tk+l f s r fS2 

w ir {sr) I ^v^Sr-i) • • • / w ix (sx)dsi •■■ dSr-tdSr (35d) 

tk Jtk Jtk 



We can restrict to % > 1 in (35a). In (35b) we can restrict to i > j + 1 as explained in 



previous subsection. In (35c), we can simplify to 

»t k 



■k+l 



v k (t)vj(t)vi(t) dt 



fc+i 



w k (t)wj(t)wi(t) dt; k>0,j<i 



Note that for the first two equalities above we need m + C(m + 1, 2) equations, where 
C(n,m) = n\/m\(n — m)\, which in total gives (m/2)(m 2 + 4m + 7). For the third one, 
we need additional m + 3C(m + 2,3). In general, it is not easy to see the formula for the 
number of equations. The number of parameters and the required degree for 0(h 4 ) are 
given by (m/2)(m 2 + 4m + 7) and iV = [(1/2) (m 2 + 4m + 5)~| respectively. 
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5 Improvements and Generalizations 



In this section, we consider techniques for improving the estimates obtained, and for 
generalizing the methods to differential inclusions with constraints. 



5.1 Improved approximate solution sets 

The previous error estimates were based on bounding the parameters appearing in the 
form of the input w(t). For example, supposing a single input v(t) G [—1, +1] and taking 
w(t) = a + a x (t - t k+1/2 )/h k satisfying J* k+1 v(t) - w(t) dt = J^ +1 tv(t) - w(t) dt = 0, we 
find |a | < 1 and |ai| < 3. However, if a = ±1, then v(t) = ±1 on [t k ,tk+i], so a i = 0- 
Similarly, if |ai| = 3 then ao = 0. 

For a given ao, we can maximise a x by taking 



v(t) 



-1 for t k < t < t k + ah k , 
+ 1 for t k + ah k <t<t k + h k = t k+1 . 



where a = (1 — a )/2. For this v, we find 



12 f tk+1 , s , s 12 



a i = J^J (t-t k+1/2 )v(t)dt = I (t-h k /2)dt- I (t-h k /2)dt 

= 3(1- (1- 2a) 2 ) = 3(1- al) 



hk i-ah 



yielding the constraint 

a 2 + \ ai \/3 < 1. 

We can therefore set 

w k {t) = a + 3(1 - al)b x (t - t k+1 )/h k with a ,h e [-!,+!]. (36) 



This will yield sharper estimates than (27). 



5.2 Differential inclusions with constraints 

Up to now, we have considered affine differential inclusions of the form 

in 

x(t) = f(x(t)) + ^g i {x{t))v i {t) with Vi G [-Vu+Vil. 

i=l 

In other words, the disturbances {y 1; . . . , v m ) lie in a coordinate-aligned box [— Vi, +Vi] x 
• • • x [— Vfc, +V^]. In many problems, the set V containing (vi, . . . ,v m ) will not be box, 
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but some more complicated set. We could use our method directly to compute over- 
approximations to the solution set by taking an over-approximating bounding box V to 
V, but this will typically yield extra solutions, even in the limit of small step size. Instead, 
we seek to restrict solutions to those of the original system. 

The right-hand-side of the differential inclusion is convex if, and only if, V is a convex 
set, so it suffices to restrict to this case. We can write 

V = {(v 1 ,...,v m )\v i e [—Vi, +VH A c(ui, . . . , v k ) < 0} 

where c : M m — > R is a convex function. (More generally, we could consider the disjunction 
of several such constraints.) The constraint c yields restrictions on the form of the u>j. 
For second-order estimates using 

I ptk+i 

w k ,i(t) = a k ,i = 7- Vi(t) dt 
rik J tk 

we simply need to introducte the constraints 

c(a k ,i, ■ ■ ■ ,a k ,m) < (37) 

at every step. For higher-order estimates, the relationship between the parameters and 
the constraint function may be more complicated; in particular, it need not be the case 
that c(w ky i (t) , . . . , w k)m (f)) < holds. 

5.3 Pseudo-affine inputs 

In this section, we consider differential inclusions of the form 

x(t)=g(x(t)) + G(x(t))q(v(t)), x(0) = x , v(t) E V (38) 

where V is compact, convex subset of W rn , and g : W 1 ->■ W n , G : W n ->• R nxp , and 
q : W m — > W. The inclusion above can be viewed in two different ways. 

One way is to consider the right-hand side as a function which is non-linear in the 
input. For example, consider a one-dimensional polynomial system with inputs, 

±(t) = X 7 vf + X vj + X 3 V!V 2 + x 5 , (v!,v 2 ) E V C R 2 . 

This has a form g(x) + G(x)q(v) by taking g(x) = x 5 , G(x) = (x 7 , x, x 3 ), and q(v) = 
(qi(v),q 2 (v),q 3 (v)) = (v\, v%, v 1 v 2 ). 

The other way is to consider the right-hand side as a function which is linear in the 
input 

x(t) E g(x(t)) + G(x(t))r(t), r(t) E q{V) 
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This corresponds to the case where, in general V, i.e., q(V) above, is convex, but not 
necessary a box as it was assumed in the previous section. For example, we can consider V 
given by constraints such as V = {v(t) \ c(x(t),v(t)) < 0} or V = {v(t) \ e(x(t),v(t)) = 0} 
for some continuous functions c(-) and e(-). 



In order to compute reachable sets of the system (38), we proceed as in the previous 
section. First we construct an "approximate" system 

y(t)=g{y(t))+G{y(t))w(t), 
and then get an error on the approximation. The local error will be essentially obtained in 



the same way as before, i.e., Theorems p 17, but with certain additional assumptions. To 



see what the assumptions should be, suppose that we want to get an error as in Theorem 



14 Then w{t) = (wi(t), ...,w m (t)) is affine and satisfies the integral equalities 



As before, we get 



q{v{t)) -w{t)dt = 

k 

tk+1 

t(q(v(t)) - w{t))dt = 0. 

tk 



I rtk+i 

a = (ai, a m ) = - J q(v(t))dt 
(6i, b m ) = — / q(v(t))(t - t k+1/2 )dt 



h 3 k 

Obviously, we can take box over-approximations for a and b, and obtain over approxi- 
mations of the reachable sets. However, if q is nonlinear, or V is not a box, but some 
general convex set, then box over-approximations for a and b could result in large over- 
approximation of the reachable sets. Therefore, if the set q(V) satisfies additonal assump- 
tions, we can get optimal results for the parameters a and b. For example, if q{V) is a 
convex set, centered around the origin, we get a 6 q(V) and b E (3/h)q(V), which gives 
optimal bounds for the coefficients a and b. 

6 Numerical Results 

We now illustrate the use of our algorithm by computing reachable sets for some simple 
systems. 
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6.1 Van Der Pol Oscillator 

We consider perturbed Van der Pol oscillator given by 



x = y 

y = —x + 2(1 — x 2 ) y + v, 



where v represents additive noise. We use the method described in Section 4.2.3 and the 



error bound (32) for additive inputs. If we take D = [0,2] x [—1,3] to be the region of 
computation, then we get K = 20, L = 31, A = 27, and H = 12. In addition, if we 
assume that v(-) G [—0.08,0.08], i.e., A = 0.08, we obtain 



27 h 

3 i ico i7l3 



\x(t k+1 ) - y(t k+l )\\ < 11.24^ + 168.17/1 



21h 



We use the algorithm described in Section |3.2| to compute the solution set for the set of 
initial points X = [0.1,0.105] x [1.5,1.505] over the time interval [0,1.5]. Because the 
bounds K, L, A, and H are rather large, we use fairly small step size, h = 0.001, yielding 
an analytical single-step error of e = 1.817092608 x 10 -7 . In Figures [l] and [2] we show 
solution set of the perturbed Van der Pol oscillator using the above values. In figure [TJ 
splitting of the domain was performed at t\ = 0.6 and £2 = 1-2. At t\ the set was divided 
in half along rr-axis, and at t% the set was divided in half along y-axis. The computed 
reachable set after T = 1.5 is a union of the following four sets: 

R(X ,T) C [1.46104,1.66704] x [-0.482307,-0.272922] 

U [1.60834,1.80823] x [-0.438931,-0.263936] 

U [1.50247,1.70832] x [-0.466819,-0.269152] 

U [1.65202,1.8518] x [-0.424135,-0.259941]. 

Moreover, if there was no splitting performed the reachable set at T = 1.5 is then 

R(X ,T) C [1.43018,1.88571] x [-0.513789,-0.197579], 

and the computed solution set is presented in [3j From the results obtained, it turns out 
that the reachable set was smaller when splitting was performed. 

Note that the set D in this case was chosen approximately, so that for initial condition 
X and time of computation T = 1.5, the solution set of the differential inclusion stays 
inside D. This is done so that analytical error does not have to be recomputed at each 
time step. In general, it is not necessary to know a-priori the region of computation. In 
fact, at each time step, we can check whether the reachable set is inside D, if not, we can 
choose new D and recompute the error accordingly. 
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Figure 1: Evolution of the Perturbed Figure 2: Evolution of the Perturbed 

Van der Pol Oscillator: splitting per- Van der Pol Oscillator: no splitting per- 

formed at t\ = 0.6 and t 2 = 1.2. formed. 

Figures [I] and [2] show that our method is effective in practice for computing rigorous 
over-approximations of the solution sets of nonlinear differential inclusions. To prove this, 
we compare the results of computation of the algorithm presented here with the ones given 
in HU. 



6.2 Perturbed Harmonic Oscillator 

The equations for the perturbed harmonic oscillator are given by 

x = y + V\ 
y = -x + v 2 , 

where v^s represent bounded noise. Suppose that the range of v± and v 2 is [—Ai, Ai] and 



[—A 2 , A 2 ] respectively. Notice that noise is additive, and therefore we can use formula (32) 
to compute the (analytical) error. In terms of our general set up we have f(x, y) = (y, —x), 
Qi = 1, for i = 1,2. Hence, we get A = 1, L = 1, H = 0, and K' = Ai + A 2 . The one step 
time error is then given by the following formula 

maxj^!, A 2 ). 



4(2 - h) h 

For comparison purposes, Table [2] is equivalent to a table given in [17] . The total time 
of computation is T = 27r, A\ = 0, and initial condition is the box (1, 0) + [—5, 5} 2 . Note 
that diameter of the set [ai, a 2 ] x [b\, b 2 ] e M 2 is max{a2 — «i, b 2 — b±}, and radius of the 
set is half of diameter. 
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From Table [2j one can see that in most cases our results are better then one obtained 
in [IT]. In case [TJ the time step is h = 2n/9 = 0.698131, for which the analytical error 
e = 0.066170 is too large to hope for sharp results. In case[3j handling the large number of 
time steps requires more sophisticated techniques for simplifying the representation of the 
intermediate sets than are currently used in our code, and this is the major contribution 
to the error. 



Table 2: Perturbed Harmonic Oscillator T = 2ix 



case 


A 2 


5 


num. of steps 


Our Diameter 


Diameter in |17j 


1 


0.1 


0.01 


9 


3.91258 


1.178825 


2 


0.1 


0.01 


100 


0.8382630 


0.8453958 


3 


0.1 


0.01 


1000 


65.4376 


0.8225159 


4 


0.1 





100 


0.8186080 


0.8253958 


5 


0.1 


0.01 


100 


0.8382630 


0.8453958 


6 


0.1 


0.1 


100 


1.018708 


1.025396 


7 


0.01 


0.01 


100 


0.1018380 


0.1025396 


8 


0.1 


0.01 


100 


0.8382630 


0.8453958 


9 


1 


0.01 


100 


8.205280 


8.273958 



When both Ax and A 2 are nonzero, i.e. A\ = A 2 = 0.1, our results and results from 
[T7] are given in Table [3j Here, we present results only for smaller time steps, even though 
in jTTj the results were given for time steps up to h = 0.799. We give both second-order 
and third-order local error estimates. We can see from Table [3] that for h = 0.25 we 
are starting to get significantly worse results then in [17J, but for smaller time steps the 
results are comparable. Here, the total time of computation is T = h (one time step), 
and 5 = 0. 



Table 3: Perturbed Harmonic Oscillator T = h 



case 


h 


Our Radius (2) 


Our Radius (3) 


Radius in [17] 


1 


0.25 


0.0420586 


0.0313667 


0.0284025 


2 


0.1 


0.0125864 


0.0108419 


0.0105171 


3 


0.01 


0.00102509 


0.00100759 


0.00100502 


4 


0.001 


0.00010026 


0.00010009 


0.00010005 



We see that the radius of the enclosure is dominated by the growth due to the noise 
in the differential inclusion. The reason why our third-order error estimates give worse 
enclosures than those of [17J is unclear; however we note that the error estimates obtained 
there were computed exactly by hand, and our automated methods are better than those 
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of [T7] based on the logarithmic norm. Moreover, in [T7] they use the 2-norm for the 
logarithmic norm which gives better results for this example. 

6.3 Rossler Equations 

The Rossler equations are given by 

= —(y + z ) + v i 

y = x + 0.2y + t> 2 

i = 0.2 + z(x — a) + v 3 

We aim to estimate the image of the initial set 

X = {0} x [-10.3 x 10~ 4 ,+10.3 x 10~ 4 ] x [-0.03 x 10" 4 , +0.03 x 10" 4 ] 

under the return map P to the Poincare section X = {x = 0, x > 0} for the parameter 
value a = 5.7 and noise Vi G [— 10~ 4 , 10~ 4 ] for i = 1,2, 3. Rather than compute the crossing 
time for each trajectory, we computed a time interval T containing the first crossing time 
by comparing the sign of x over the sets Rk, and used the estimate {0} x P(Xq) C 
R(X ,T). 

With time step h = 0.005, total time T = 11.1, and region of computation D = 
([-25,25], [-25,25], [-25,35]), we obtain an analytical error of e = 8.586 ■ 10" 8 and 

R(X , T) = ([-0.15572, 0.15391], [-3.75926, -3.41772], [0.03139, 0.03398]). 

In [T7], R(X ,T) = ([-0.211150,0.20888], [-3.69781,-3.47352], [0.03117,0.03327]). (They 
did not specify the time step or the total time it took to compute the value of the poincare 
map R(Xq,T).) In this case neither of the sets is better then other, but they are com- 
parable, and hence we show that our algorithm can also provide good estimates when 
computing over rather difficult regions, see [T7] . 

7 Concluding Remarks 

In this paper, we have given a numerical method for computing rigorous over- approximations 
of the reachable sets of differential inclusions. The method gives high-order error bounds 
for single-step approximations, which is an improvement of the first-order methods pre- 
viously available. By providing improved control of local errors, the method allows for 
accurate computation of reachable sets over longer time intervals. 

We give several theorems for obtaining local errors of different orders. It is easy to 
see that higher order errors (improved accuracy) require approximations that have larger 
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number of parameters (reduced efficiency) . The growth of the number of parameters is an 
issue, in general. Sophisticated methods for handling this are at least as important as the 
single-step method. The question remains as to approximate solution (Theorems [8 17) 
yields the best trade-off between local accuracy and efficiency for computing reachable 
sets. The answer is not straightforward and most likely depends on the system itself. 
In future work, we plan to investigate the efficiency of the algorithm on the number of 
parameters for various examples. 

We have only considered differential inclusions in the form of input-affine systems, and 
give a brief sketch of how these methods can be applied to other classes of system. We 
also plan to provide a more detailed exposition of the method in these cases. Moreover, 
the local error that we obtain is a uniform bound for the error in all components. It 
should be possible to give slightly better componentwise bounds. 
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